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In this work we study the dynamics of free 3D relativistic Gaussian wave packets with different 
spin polarizations. We analyze the connection between the symmetry of initial state and the dy- 
namical characteristics of moving particle. The corresponding solutions of Dirac equation having 
different types of symmetry were evaluated analytically and numerically and after that the electron 
probability densities, as well as, the spin densities were visualized. The average values of veloc- 
ity of the packet center and the average spin were calculated analytically, and the parameters of 
transient Zitterbewegung in different directions were obtained. These results can be useful for the 
interpretation of future experiments with trapped ions. 
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I. INTRODUCTION 

The Dirac equation belongs to the most important 
equations in modern physics. It predicts the existence 
of electron spin and magnetic moment, gives the nat- 
ural description of the positron states, describes fine 
structure of the energy spectrum of hydrogen-like atoms. 
The quantized solutions of Dirac equation are consid- 
ered to be a natural transition to quantum field theory. 
Moreover, the one-particle relativistic quantum mechan- 
ics describes the unexpected electron dynamics including 
Schrodinger's Zitterbewegung (ZB)i and Klein paradox.^ 
The trembling motion of relativistic particles or ZB is 
caused by the interference between positive and negative 
energy states which form the electron wave packet. The 
frequency of ZB is determined by the gap between two en- 
ergy bands, and the amplitude of oscillation of the wave 
packet center is of the order of the Compton wavelength. 

The results of the first experimental observation of ZB 
phenomena were published recently in the paper by Ger- 
risma et. al^ For the ZB simulation the experimental- 
ists used a linear Paul trap where ion motion can be 
described by one-dimensional Dirac equation<^ The au- 
thors of Ref. [3] study the motion of Ca~^ ion and de- 
termined its position as a function of time for different 
initial conditions. As was shown in Ref. [4] the solution 
of the 3-1-1 Dirac equation can also be simulated using a 
single trapped ion with four ionic internal states. In this 
case the ion position and momentum are associated with 
respective characteristics of 3D Dirac particle. 

The dynamics of relativistic one-dimensional wave 
packet for the first time was investigated numerically by 
Thalleri^ He visualized the initially localized solutions 
of single particle Dirac equation and observed the trem- 
bling motion of the wave packet centers as well as some 
other phenomena which are caused by the interference 
of positive- and negative-energy states. J. Lock found 
that ZB oscillations of localized initial states have tran- 
sient rather than sustained character.^ Some interesting 
examples of the relativistic dynamics of 3D electron wave 
packets were presented in Ref. [7]. It should be stressed 



that ZB oscillations exist only in the one-particle rela- 
tivistic theory. Krekora, Su, and Grobe demonstrated 
analytically and numericallj*^ that quantum field theory 
does not permit any Zitterbewegung of real electrons. 

The oscillatory ZB motion of 3D electron wave packets 
in crystalline solids for the first time was predicted in 
Ref. [9]. This phenomenon has been considered for 2D 
electron gas with Rashba spin-orbit coupling in Ref. [10, 
11], in narrow gap semiconductors in Ref. [12], in carbon 
nanotubesi^, in single and bilayer graphenei^"— and also 
in superconductorsi^. 

In the present work we investigate the relativistic dy- 
namics of 3D wave packets. In Sec. II the general prop- 
erties of symmetry of Dirac equation solutions which de- 
termine the dynamics of wave packet are analyzed. After 
that, in Sec. Ill we consider the evolution of initial states 
with different initial symmetry. The electron probability 
densities evaluated analytically and numerically are visu- 
alized. The directions of average velocity of wave packets 
as well as the phenomenon of high-frequency Zitterbewe- 
gung for different initial polarization are determined. In 
Sec IV the symmetry and structure of spin densities for 
relativistic packets are analyzed. Finally, in Sec. V, we 
conclude with the discussion of results. Some auxiliary 
results are found in Appendixes A and B. 



II. THE SYMMETRY OF THE RELATIVISTIC 
WAVE PACKETS 

In this section we shall first study some general sym- 
metry properties of the electron wave function which de- 
termine the dynamics of wave packets. So, we start with 
the famous Dirac equation for a four-component wave 
function \1> = (^i, \i;'2, ^3, ^4)^ 



where 



5* 

ih— = i/*, 

dt 



H = —ic:haW + fimc , 



(1) 
(2) 



c is the light velocity, m is the electron mass and 
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(3) 



The four independent free-particle solutions of Eq.(l) 
for given momentum p and energy _E, can be written in 
the form 



M/,>(f,i)=c-''^*/Vp(?)Ur(p), 1 = 1,2,3,4, (4) 



where ipp{f) = j^:;^^e'P'/^\ E = ±Ajr, \f; 
\J'P^c?' + vrP-d^ , and Ur{p) are the free Dirac spinon 
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Ui{p) = N 




U2{p) = N 




-P3l 







E = -Xf;<0, 
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Ap- + mc^ 



Here we use the normalization condition 

U+{p)Ur'ip)^Srr', 

which means 



N^ y'(Aj? + mc2 )/2Ap-. 



(7) 



(8) 



(9) 



Let us consider now the symmetry of the solutions of 
Dirac equation with respect to the space reflection. Sup- 
pose that initial wave packet is symmetric: |\l/(r, 0)p = 
\'^{—r, 0)p. As it known this symmetry is conserved, i.e. 
|vl>(f^i)|2 = |\]>(_f^i)|2 Qj^jy if the full function including 
the spinor part is the eigenfunction of the parity operator 
P — /3R, where the matrix /3 is determined by Eq.(3) and 
R is the inversion operator for f: R'^{r,t) = ^{—r,t). 
Also the Dirac equation is invariant under the replace- 
ment 



is the corresponding component of spin operator E 



a 



So, if the wave function satisfies (in particular 



at i = 0) the relation 

then the parity in x, y-plane is conserved. 



(11) 



(12) 



Besides, we may determine the reflection transform z — >■ 
—z in Hilbert space as the operator P^ = YizPRz which, 
just as the operators P and Px,y, commutes with Dirac 
Hamiltonian (2). If the initial wave packet has the certain 
party z — > — z 



P,*(f,0) = T*(r,0), 



(13) 



then ^(r, t) satisfies this equation for all times that leads 
to the conservation of the symmetry 



\^{x,y,z,t)\^ = \^{x,y,-z,t)f 



(14) 



Similarly we may introduce the operators P^ = Y,xPRx 
and Py = Y,yj3Ry connected with the reflection trans- 
forms X — > —x and y — > —y correspondingly. 

Thus the symmetry of the solution ^(r, t) with respect 
to full (r — > —r) or partial space symmetry depends not 
only on the space symmetry of the initial wave function, 
but also on the ratio between its components. This state- 
ment is valid with respect to other types of symmetry. 
Below we are interested in the dynamics of the initial 
wave packet of the form 



*(f,0) = 



F{r) 



E Wi 




(15) 



where ipi are the complex numbers and the space part 
F[r) satisfies the normalization condition 



I \F{r)\^dr = 1. 



(16) 



Specifically we suppose that the probability density at 
t — |^(r, 0)p = |F(r)p is spherically or axially sym- 
metric. So let F{r) has the form 



*(x,2/,z,i) ^ *'(r,i) = Px,y^{r,t), 



(10) 



where Px.y = Y^zRxRyi Rx{Ry) is the inversion operator 
for x{y)- component, Rxf{x,y,z) = f{—x,y,z) and S^ 



F{r)^F{p,z)e" 



(17) 



where p, z, a are cylindrical coordinates and m is an 
integer. The considered wave packet remains an axial 



symmetrical (in x, y-plane) for all the times if its spinor 
part (p is one of the eigenfunctions of operator E^ 



a first example let us consider the time evolution of the 
initially localized wave function of the form 



^ ^ {(pi (p3 0)^ ^^ ^ ^(0 ^2 ipif 



VI'^iP + l'^sl 



V\f2? + \f4\' 



Indeed it is readily to show that the function 



*(f,0)=i^(/9,z)e™"$i, 



(18) 



(19) 



is the eigenstatc of z-component of total angular momen- 
tum operator I^ — Iz + ■^^z 



*(f,0) 



V2 



(24) 



As was shown above (see Eq.(18)) for such polarization 
of the wave packet the probability density conserves its 
axial symmetry if the space part of the wave function 
(24) F{r) has the form (17). 

The wave function ^(r, i) can be expanded in plane 
waves 



/^*(f,0) = (TO+-)*(r,0), 



(20) 



where Iz — —id/da. 

For time i > the general expression for the wave 
function is 



*(f,i)= / Lpp{r)^{p,t)dp. 



(25) 



In the momentum space the components of the bispinor 
wave function 



*(r, t) = (*i(p, a, z, t), *2(p, a, z, t), 

*3(p,a,2,i),*4(/5,a,z,i))'^. (21) 

But Iz is a conserved quantity so that the ^(r, i) obeys 
Eq.(20) too. Solving this equation we find the a - 
dependence of components ^i(p, a, z, t) 

*i = e™"fi(p,z,t), *2 = e'('"+i)"f2(p,z,t) 

*3 = c'"^"f3(p,z,t), vl/4=c'(-+i)"f4(p,z,t), (22) 

which immediately leads to conclusion that the probabil- 
ity density is axially symmetric in x, y-plane. 



\^ir,tr==J2\MP,^,t)? 



(23) 



Note that the initial state *(r, 0) = F{p, z)e™"$_i is 
the eigenfunction of Iz too with an eigenvalue m ~ 1/2. 
Notice also that as would be shown below the spherical 
symmetry of the initial wave packet, Eq. (19) (with m — 
and F{p,z) = F{r), r — ^ p^ + z^) turns into axial 
one. The above statements concerning the wave packet 
symmetry arc justified by our analytical and numerical 
calculations for different initial conditions. 



III. DIRAC WAVE PACKETS DYNAMICS 



*(p, t) = (*i(p, i), *2(P, t), *3(P, t), *4(P, t)f, (26) 

are given as a linear superposition of the positive- and 
negative-energy solutions (5) and (6) 



vI/(p,f) = ^a(p)t/.(p)e-'^''/^ 



(27) 



r=\ 



where _Bi_2 ~ Xp and E34 ~ —Xp- Here C'r{p) is to be 
determined by the Fourier expansion of ^(r, 0). Straight- 
forward calculation using (5), (6) and (24) gives 



*i(p,i) 



m 

2V2 



-.-i^pt/fin ^ 



mc + cp3 , 



h 



_^giApt/h(-j^ 



*2(p,i) = *4(p,t) = - 



mc" + cp3 , 



*3(p,i) 



m 

2V2 



g-iAjjt/H^2 — 



mc — cp3 , 



Ap 



+c'^p*/'^(l + 



mc — cp3 , 



(28) 



^:2f(MPl±Msu.hEL, (29) 
Ap-\/2 ri 



(30) 



i) The propagation of cylindrically symmetric 
wave packet 

Now we describe some peculiarities of the striking kine- 
matics of three-dimensional relativistic wave packets. As 



Here f{p) is the Fourier transform of the function F{r). 
To obtain the components of wave function ^^(r, f) we 
substitute Eqs. (28)-(30) into Eq. (25). Then for the 
axially symmetric initial wave packet F(f) = F(p,z), 



where p"^ = x'^ + y^ , i.e. /(p) = f{pj_,Pz), p\ = pI+pI^ 
after integrating over the angular variable we shall have 



*i(p,z,t) 



1 



+ CX3 



4V7i^ 



e K dpz / f(p±,pz)p±x 



iApt/H(l + "1^ +^P3 ) 



+e*Ap*A(l-!^^^±^) 



Ap* 



Ap 



M^)dp^, (31) 



+ 00 



*2(p, a, 2^, i) = *4(p, a, z, t) = 



e K dpzX 



2V^ 

— OO 

x|/(p^,p.)^Ji(^)sinMdp^, (32) 



]1|Ii:t,y.!aCtiarj)|^ 




^3(P,^,0 = 



1 



+ 00 



— ex 



e K dpz / f(p±,Pz)pj.x 





mc — cp3 



Ap 



)+ 



+e'^p*/^il + 



m(f - cp3 , 



M^)dp^, (33) 



where Jo(u), Ji{u) are Bessel functions. Using these ex- 
pressions one may find the dependence of the electron 

4 

probability density ^ |^i(^, OP on coordinates and time 

i=l 

for the wave packet (24) . To illustrate our results in this 
Section and everywhere below we choose the function in 
the Gaussian form 



F{r) = 



1 



dVATT^ 



= e 2d'= 2A- 



■+ikoz 



(34) 



where d and A determine the width of wave packet in 
the x, y-plane and in z direction correspondingly, and 
Pz — fi-ko is an average momentum of the wave packet. 

In Fig. 1(a) we plot the electron probability density at 
z = for Gaussian wave packet with initial momentum 
ko = and A = 5, d = 1 at time t = 7.5. Here and 
below all distances are measured in units of the Compton 
wavelength Afe = h/mc, the time is in units of to — Xk/c. 
This result was obtained by using the numerical method 
("leap-frog" algorithm) described in Appendix A. Just as 
we expected Fig. 1(a) demonstrates the axial symmetry of 
the considered distribution in a;, j/-plane; the probability 
density has a form of cylindrical wave propagating from 
the point of origin and having some maxima. 



FIG. 1: (Color online). The electron probability density for 
the initial Gaussian packet, Eqs. (24), (34): (a) at 2; = 0, 
fco = and A = 5, d = 1 at time t = 7.5; (b) at y=0, fco = 
and A = 5, d = 5 at time t = 7.5. 



To analyze the motion of the packet we have to find 
the average value of velocity of the packet center. In the 
momentum representation 



V^{t)=c hi'+{p,t)a,^{p,t)dp, i = 1,2,3. (35) 

Substituting Eqs.(28)-(30) into Eq.(35) and using the ex- 
pressions for matrices a^ (Eq.(3)) we obtain 



Mt) = c dp\f{p)\ 



2— (1-cos ^)4-;- sm — ^ 






Ap 



n 

(36) 



Vy(t)^c dT^f{p)\ 



(?P2Pi,^ 2Ap-i cpi . 2\pt 



n 

(37)- 



Vz{t)=c dp\f{p)\ 



''^''(l-^)cos 



Z\ftt 



XI 



XI 



, (38) 



Since the Fourier transform f{p) of Gaussian wave packet 
(34) has an axial symmetry in px,Py- plane 



d^/d 



(39) 



the components of velocity T4 = Vj, = as it follows from 
Eqs.(36),(38). Otherwise, owing to the axial symmetry of 
spacial distribution of the electron density, the average 
coordinates of packet x = y = 0. As a result, mean 
components of velocity in x, y-plane are equal to zero. 

The motion of the packet center in z-direction expe- 
riences rapid oscillations commonly known as Zitterbe- 
wegung (the second term in square brackets in Eq.(38)). 
Besides, the wave packet displaces slowly with constant 
velocity Vzo (the first term in square brackets in Eq.(38); 
see also Eq.(B.8)) even if its momentum p^ = hko = 0. 

The existence of the constant component of average ve- 
locity Vzo in this case can be understood from the relation 
between velocity and momentum depending on the sign 
of energy: for the wave packet consisting of the states 
with positive (negative) energy a positive momentum pz 
corresponds to a positive (negative) velocity Vz- Let us 
represent the time- independent probability density in the 
momentum space |^(p)p as a superposition of the pos- 
itive energy part ^+(p) and the negative energy part 
*-(p) 




FIG. 2: (Color online). The dependencies W±{pz) of positive 
(solid line) and negative (dashed line) energy parts for Gaus- 
sian initial wave packet, Eqs. (24), (39) for fco = 0, A = 5, 
d — 1 (a), and for fco = 1, A = 5, d = 5 (b). 



\^ip)\' = \^+{p)\' + \^^{p)\' 



where \^±ip)\^ = EI*tt(P,OI', *.±(P,t) « 
Using Eqs. (28)- (30) one may readily show that 



(40) 



aTl- 



\^±{P}? 



urn' 



(i± 



CPz 



(41) 



Fig. 2 shows the dependence W±{pz) 



2t: J \'i'±{p)\~pj_dp± (in arbitrary units) on pz (in 



units of mc) for the Gaussian initial wave packet, Eqs. 
(24), (39) for ko = (a) and for ko — 1 (in units of 
mc/h) (b). 

We see that for the case fco = the momentum dis- 
tribution for positive- and negative- energy components 
is shifted towards positive and negative momentum pz, 
respectively. But in both cases such dependence leads to 
the motion of each parts (and the whole wave packet) 
with a positive velocity along z-axis. We also see in 
Fig. 2(a) that the function W+{pz) has essential overlap 
with the function W-{pz) that is a necessary condition 
for existence of Zitterbewegung of the packet center (see 
Eq.(38)). The value of constant component of velocity 
Vzo depends on the ratio between the initial width of the 
wave packet and the Compton wavelength. In the limit- 
ing case d >> 1 t4o is much less than the light velocity: 
VzQ — f (^)^- In particular, for the symmetrical wave 
packet (fco = 0) of width d = A = 5 t^o « 0.02c. 

Let now the initial wave function F{r) describes the 
wave packet moving along z-axis with average momen- 
tum Pz = hkQ. Then the distribution of the full electron 
density in x, y-plane is similar to one shown in Fig. 1(a). 



However the essential difference appears in the charac- 
ter of evolution of the wave packet in x, z (or y, z)-plane 
(Fig. 1(b)). The dependencies W±{pz) for the states with 
positive- and negative- energy parts for this case is shown 
in Fig. 2(b). We see that both components consist of 
positive momentum pz . For the smaller negative-energy 
parts, this corresponds to the motion with negative ve- 
locity along z-axis. Thus in the position space the initial 
wave packet splits into two packets propagating in the 
opposite directions along z axis. 

ii) Asymmetrical wave packet evolution 
We next consider another example of the initial spin 
polarization of the electron 



*(f,0) 




(42) 



where as before the function F{'r) is determined by 
Eq.(34). 

Note that the example under review is invariant with 
respect to the reflection transformation z — >■ — z, i.e. the 
expression (42) satisfies Eq.(13). As was shown above 
this means that the probability density is an even func- 
tion of z at all times. Performing the same kind of calcu- 
lations as for symmetrical wave packet we find the com- 
ponents of the initial wave function (42) at i > in the 
momentum space. 



*i(p,0 



m 

2V2 



-iXpt/h/-^ I 



+ c{pi -ip2) 
Ap 



)+ 



_^eiApt/H(^-^ 



mc^ + c(pi - ip2) 



(43) 



*2(p,t) = -*3(p,<) = 



icf{p)P3 . >^pt .... 



*4(p,i) = 



2^2 



-iApt/H^2 



mc^ -c(pi +ip2) 



iAgt/H^i , mc"-c(pi+ip2) 



+e'^p*/"(l + 



Ap 



)+ 
(45) 



The components of *(r, t) can be obtained directly by 
the Fourier transform of Eqs.(43)-(45). 



*i(p,a,i) 



1 



+ C» 



2V7rr 



e K dpz / f(p^,pz)p±x 










-— — sm^Jo(^) 



dp±, 



(46) 



+ 00 



*2(/5,Z,t) = -*3(P,2:,0 



2V7r?i^ 



p^e K dpzX 



OO 

X J f{p^,p.)^^sin^M^)dp^, (47) 



lf*».y.F<;i-:™)l 




FIG. 3: (Color online). The electron probability density for 
the initial Gaussian packet, Eqs. (42), (34): (a) at z = with 
kg = and A = 5, d = 1 at time t = 2it; (b) at y=0, ko = 1 
and A = 5, d = 2.5 at time t — 7.5. 



V^t) 



dp\f{p)\' 



9 9 9 9 

^+(1-— %i)cos 



Z\i^Z 



X2 



X2 



(49) 



Vy{t)=C dp\f{p)\ 






mc 2\0t 

'— — sm — i— 



h 

(50) 



*4(p,a,i) 



1 



+ 00 



2V7rr 



e B dpz / f(p±,pz)p±x 



cos -^ Jo(-^) + -— sm -^e'"Ji(^)- 
h h Aff Ti n 



-^— sm^Jo(^) 



dpA 



(48) 



Using these expressions one may find the full electron 
density |^(r, i)p. Figure 3(a) shows the corresponding 
probability density distribution in the plane z = at 
time t — 27r. The parameters of wave packet are: kg = 
and A = 5, d = 1. Comparing Figs. 1(a) and 3(a) we 
see that the change of initial spin polarization leads to 
the fact that the wave packet being axially symmetric in 
X, y-plane at i = loses its symmetry at t > 0. 

The kinematics of the wave packet can be characterized 
by the average velocity of its center with components 



Ut) = c dp\fip)\ 



12C -P1P3 
^1 



1 — COS ■ 



ZAift 



(51) 



As in the previous case the wave packet center drifts with 
constant velocity (the first term in square brackets in 
Eq.(49)), but now it is directed along x axis. Besides, 
one can see from Eqs. (49) and (50) that the packet 
center performs damped oscillations (Fig. 4) along x and 
y directions. We also see from Eq.(51) that Vz — that 
is the result of the symmetry under the replacement z — >■ 
-z, Eq.(14). 

We now consider the behavior of the initial wave packet 
with nonzero momentum pz = fik^. Obviously such ini- 
tial state is not an eigenfunction of the parity operator 
Pz- Nevertheless, the probability density remain to be 
a symmetrical function relatively to the reflection trans- 
form z — > —z. Indeed, one may check that the Dirac 
equation (1) is invariant under the transformation 

^(r, t) -^ *'(a;, y, z, t) = a^-^*{x, y, -z, t). (52) 




FIG. 4: (Color online). The average projections of velocities 
t4/c (a) and Vyjc (b) versus time t (in units of \k/c) for 
Gaussian initial wave packet, Eqs. (42), (34) with fco = and 
A = 5, d= 1. 



So that if the initial wave function satisfies the equation 



,**(a;,2/,-z,0) = *(a;,y,z,0). 



(53) 



then, as one can check using Eqs. (46)- (48), this relation- 
ship is valid at i > and consequently Eq. (14) holds. 

The character of the motion of the packet center in 
x, y-plane is similar to the case k^ — Q. Namely, the cen- 
ter of the wave packet drifts along x-direction with con- 
stant component of velocity and oscillates along x and y 
axis (Zitterbewegung). As it known the ZB is significant 
if the subpackets with positive and negative energy have 
overlap in the position space. But if the initial average 
momentum of the wave packet is nonzero, both subpack- 
ets move with the opposite velocities along z axis that 
leads to their spatial separation. So, the amplitude of 
the ZB decreases more rapidly than for the case ko — 
(compare Fig. 4 and Fig. 5). Notice that this result 
is valid for the narrow enough wave packets with width 
d < 1 and dko ^ 1 (for very large d the ZB oscilla- 
tions are almost undamped). The constant component 
of packet center velocity V^o also depends on the initial 
momentum hko. In fact, as it follows from Eq. (49) it 
decreases as fco increases. 



IV. SPIN DYNAMICS 

At present we consider the average value of the spin 
operator E for Dirac particle 



4 

S^(i)= I dr^^{r,t)^ f drJ2'^+ir,t)j:^^,ir,t), 

(54a) 




I 10 n 14 



FIG. 5: (Color online). The average projection of velocity 
t4/c(a) versus time t (in units of \k/c) for Gaussian initial 
wave packet, Eqs. (42), (34) with fco = 1 and A = 5, d = 2.5. 



or in the momentum representation 

t^{t) = J dpY, *+(p, i)s^*.(p, t). 



(546) 



One can verify that the spin densities S^(r, t) for the 
arbitrary wave function with components \l/i = 4'i(r, i) 
are 



E,(f,i) = 2Re(*J^2 + *3*4), 



Sj,(f,t) = 2Im(*J^2 + *;^4), 



E,(r,t) = |vI/i|2-|M/2| 



i*.^r-i*4 



(55) 



(56) 



(57) 



i) Cylindrically symmetric wave packet 

The form of Eqs.(55)-(57) and the expressions (31)- 
(33) for the components of wave function show that for an 
axially symmetric wave packet, Eq.(24) the component 
of spin density Sz is an axially symmetric too both for 
fco = and fco ^ (Fig. 6(c)). One can also verify that 
Y,x and Ej, components can be written in the form 



X- g{p,z,t), 



y-9{p,z,t), (58) 



where the expression for the function g(p, z, t) is enough 
cumbrous and will not be presented here. Thus the x, y- 
plane projection of vector S(r, t) is directed along the 
unit vector of cylindrical system Cp. The spin densi- 
ties T,x{x, y, 0, t) and Ey(a;, y, 0, t) are represented in Fig. 
6(a),(b) for the Gaussian wave packet with the average 
momentum pz = ?ifco = 0. We see that 'E,.j;{x,y,Q,t) 
(Sy (x, y, 0, t)) is the antisymmetrical function oi x{y) and 
the spin density Sz(a;,y,0,t) conserves its axial symme- 
try. 

The average values of spin operators for this polariza- 
tion can be found by using Eqs. (28)- (30) for wave func- 
tion in the momentum representation and previous defi- 
nition, Eq.(54b). 
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FIG. 7: (Color online). The distributions of the components 
of spin density Ea;(a;,0, 2, t), 'Sy{x,0,z,t) for initial Gaussian 



packet, Eqs. (42), (34) with fco 



and A = d = 3.64 at time 



FIG. 6: (Color online). The distributions of the components 

of spin density E^(a;,y,0,t),Sa(a::,y,0,t), S^x, J/, O,i)forini- that together with Eqs.(55)-(57) leads to the result 

tial Gaussian packet, Eqs. (24), (34) with fco = and A = 5, ^.^(^x, y, 0, t) = T.y{x, y, 0, t) = E^(x, y, 0, t) = 0. 

d = 1 at time t = 7.5. rpj^g analysis of expressions (46)-(48) and (57) shows 

that in x, z-plane the 2;-coniponent of spin density is an 
antisymmetric function of z that is connected with dis- 
cussed above symmetry of probability density with re- 
spect to the replacement z ^- —z. 
dp, Fig. 7 illustrates the distributions of Y,xix,0, z,t) and 
T,z{x,0, z,t) for Gaussian wave packet with fco — 1 a-nd 
A = d = 3.64, at i = 8. 

By inserting Eqs.(43)-(45) into Eq.(54b) we find for 
the average components of spin 



s,(t) 



ty{t) = 






cp2>yp 



sm -l-c p3Pi(l— cos ) 



(59) 



\mi 



. ZXpt 2 

-cpiAp-sm -^ ' " 



+c P3P2(1- 



-cos ■ 



2 Apt 
(60) 



dp, 



t^t) = 



\mi 



2A-t 
mc^P3{l- cos ^-)dp, (62) 



s.W 



\mi 



[XI-^pD+c^pI^os^-^ 



dp. (61) 



ty{t) 



\l{P)? 



Xf! 



-cps sm 



^dp. 



(63) 



As it follows from these relations only T,z{t) is not equal 
to zero for the considered wave packet. The first term 
in square brackets of the last formula corresponds to 
the constant component of ^z{t) and the second one de- 
scribes the typical transient Zitterbewegung. 

ii) Asymmetrical vi^ave packet 

One can easy verify that the second example, Eq.(42) 
considered in our work corresponds to the initial spin 
density S(r, 0) = 0. Really, as it follows from Eqs. (46)- 
(48) the components of wave function for the packet with 
fco = at z = obey the relations 



*t(a;,y,0,i) = *4(a;,y,0,i), 

*;(a;,2/,0,i)=*3(a;,y,0,t)=0. 



S,(t) = 






2Ap*t mc^pi 
cp2 sm — r^- — I ^ X 



2Ap't 
x(l — cos— -^) 



h 

dp. 



Ajj 



(64) 



Obviously for a symmetric wave function f{p) — f{—p) 
(that means fco = 0) these expressions give E^; = Sj, = 
fiz = 0. Note that Sz remains to be equal to zero also 
for pz = hko ^ 0. It may be shown that in this case 
Eqs. (62), (63) describe the spin "precession" (in x, y-plane 
about vector fco) which has a transient character. Such 
phenomenon for a hole system, described by Luttingcr 
model was discussed in the recent work of authorsJ^ 



V. SUMMARY 

In this work we have studied the quantum dynamics 
of relativistic particles represented by three-dimensional 
Gaussian wave packets with different initial spin polar- 
izations, described by the Dirac equation. The analysis 
of the general symmetry properties of solutions of one- 
particle Dirac equation allows to predict the direction of 
average electron velocity as well as the direction of trem- 
bling motion. In particular, the evolution of spherically 
and cylindrically symmetric initial Gaussian wave packet 
demonstrates that the wave packet with initial polariza- 
tion, which is determined by bispinor (1 1 ) has 
cylindrical symmetry at all times i > 0, but the wave 
packet with initial polarization (1 1 ) loses its 
cylindrical symmetry at time t > 0. The influence of the 
symmetry of initial wave packet on the distribution of 
spin densities is analyzed. 



Appendix A 

The "leap-frog" algorithm^^ is applied in a spatial grid 
of bin-sizes Ax, Ay, Az and with time step At: 



*(f,t + 2At) = *(f,t)-2iAti/(f)*(r,i + At), (^.1) 

The spatial derivatives are computed symmetrically. Re- 
flecting boundary conditions are applied on a very large 
grid (running stops before reflections occur if necessary) . 
We use the norm as the stability measure of the algo- 
rithm (1). In accordance with von Neumann stability 
analysis^S (for large component plane waves) the stabil- 
ity region (d-spatial grid bin, At-time step) is: 



1 - (Atf - 2{dAtf - A{Atf > 0, 



{A.2) 



Thus, for a single precision calculation the loss of norm 
can be kept within 10^^-10~^ in a 10'^ time step run. It 
should be noted that in the case of Zitterhewegung, i.e. 
of the spatial oscillation of the wave packet, one more 
condition have to be imposed to the lattice sizes: 



Aa; -- Ay -- Az < — , 
mc 



(A.3) 



The conditions (A.2) and (A. 3) were fulfilled in all our 
calculations. 



relativistic kinematics. In particular, besides the rapid 
oscillations (ZB) the wave packet can drift with constant 
velocity although its average momentum is zero. One 
can show that in this case the direction of such motion 
coincides with the direction of initial average velocity. In 
fact, for the second example considered in this work (Eq. 
(42)) only V^{t = 0) = c and Vy{t = 0) = V^{t = 0) = 
0. Therefore such initial spin polarization leads to the 
motion of the wave packet with constant velocity along 
X axis (see Eqs. (21), (26)). It is easy to see that in the 
first example, Eq. (10), the wave packet at i = has 
the velocity directed along z axis, so the motion at t > 
occurs in this direction. 

Let now find the drift velocity of particle for the case 
of arbitrary initial polarization 




{B.l) 



El 



^{p,t)^Af{p) 



where i^i are the complex coefficients, A 



and f{'p) is to be determined from the Fourier expansion 
of coordinate wave function F{f). (We do not suppose 
that F{r) and /(p) have any symmetry). 

At t > the wave function in momentum representa- 
tion is 



^ip,t) = ^+{p,t) + ^^{p,t) = 

2 4 

Y, C,(p)t/.(p)e-'^p'/^ + J2 Ci(p)Ui(p)e'^p*/^ (B.2) 

2=1 i=3 

Using the Eqs (5), (6) and (8) for free Dirac spinors we 
find from Eq. (B.2) 



Ci = Af{p)Niipi + 7(pi - ip2)ip4: + IPsVi), 
C2 = AJ{p)N{ip2 + 7(pi + ip2)^i - IPz^i), 
C3 = Af{p)N{-^p3ipi - 7(pi - ip2)'^2 + (^3), 
Ci = Af{p)N{-j{pi +ip2)^Pi +W3V2 +'^4)-(B.3) 

The density of /i-component of velocity {fi — 1,2,3) in 
the momentum space is 



Vf,{p,t)^^+{p,t)%^{p,t), V^, 



CUu 



{BA) 



Obviously the time- independent part V^o of V^ {p, t) is 
defined as 



Appendix B. Drift velocity for the arbitrary initial 
wave function 



As was shown in previous investigations the motion 
of the Dirac wave packet center does not obey classical 



^mo(p) = *T(P' i)^M*+(P, t) + *=(P, t)V^^-{p. t), 

{B.5) 
where ^+(p, t) and ^_(p, t) corresponds to the contribu- 
tion of positive and negative energy into Eq.(B.2). One 
may check that 



10 



£PiiA.. i 7 - 1 2 



So, using Eqs. (B.2), (B.5) and (B.6) we obtain 



(B.' 



V,o= /%^'-'^ 



i\Ci\ 



IC2P 



IC3P '^'^ 



C4ndp. (B.7) 



Substituting the expression for the coefficients Ct, 
Eq.(B.3) into the Eq.(B.7), we find the constant velocity 
of wave packet center 



v;.o = 



E I'Pi 



\fip)\ 



^1 



(Note that in this equation there is no summation over 
fj,). Let the initial Gaussian wave packet be spherically 
symmetric, i.e. f{p) is determined by Eq.(39) with d — 
A. Then if the average momentum hko — the second 
term in Eq.(B.9) equals to zero and the value of integral 
in the first term does not depend on fi. Thus the direction 
of the constant velocity of wave packet coincides with 
the initial one. In the case kod >> 1 as it follows from 
Eq.(B.9) Vzo >> Vxo, Vyo and the asymptotic direction of 
the average velocity is along z-axis, i.e. along the average 
momentum pz = hko. A similar result was obtained in 
our previous worki^, concerning the propagation of the 
wave packet in graphene. 

Note that Eq.(A.8) is valid also for the most general 
initial wave packet of the form 



mc(|<^ip+|<^2p-|<^3p-|</^4n 



*(p,0) = (*i(p),*2(p1,*3(p),*4(p))' 






if we rename ^i{p} 



f{p)Vi 



in Eq.(B.8) . It is not 



difficult to check that the expressions for the constant 
components of velocity obtained earlier for the examples 



J, ■ • - i J- ^u- • • ii, Eq. (24) and Eq. (42) follow also from the general equation 

it IS convenient to represent this expression using the /ji q\ ' ^ ^ ' o -^ 



initial components of velocity V^(0), so that 



(B. 
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V,o = V,{0) I |/(p)P^dp- 



{\^i\ 



l/(p)P|x 



ml 



m\ 



m\ 



-dp. 
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